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ABSTRACT 

The injection of turbulence in the circum-galactic medium at redshift z = 2 is in- 
vestigated using the mesh-based hydrodynamic code ENZO and a subgrid-scale (SGS) 
model for unresolved turbulence. Radiative cooling and heating by a uniform Ultravi- 
olet (UV) background are included in our runs and compared with the effect of turbu- 
lence modelling. Mechanisms of gas exchange between galaxies and the surrounding 
medium, as well as metal enrichment, are not taken into account, and turbulence is here 
driven solely by structure formation (mergers and shocks). We find that turbulence, 
both at resolved and SGS scales, impacts mostly the warm-hot intergalactic medium 
(WHIM), with temperature between 10^ and 10^ K, mainly located around collapsed 
and shock heated structures, and in filaments. Typical values of the ratio of turbulent 
to thermal pressure is 0.1 in the WHIM, corresponding to a volume-weighted average 
of the SGS turbulent to thermal Doppler broadening 6t/&therm = 0.26, on length scales 
below the grid resolution of 25 kpc h~^. In the diffuse intergalactic medium (IGM), 
defined in a range of baryon overdensity S between 1 and 50, the importance of turbu- 
lence is smaller, but grows as a function of gas density, and the Doppler broadening 
ratio is fitted by the function ^t/^thcrm — 0.023 x S'^-^^. 
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1 INTRODUCTION 

The properties of the diffuse gas component, not yet ac- 
creted in collapsed cosmic structures, are subject of study 
in order to understand the census of the baryonic matter 
which d oes not be long to galaxies. Early numerical investi- 
gations (!Cen fc Ost riker 1999; Dave et al. 2001) have shown 
that a substantial fraction (30 to 40 per cent at redshift 
2: = 0) of this gas resides in a baryon phase with tempera- 
tures between 10^ and 10^ K, the Warm-Hot Intergalactic 
Medium (WHIM), and a similar a mount is contained in a 
colder (T < 10^ K) diffuse phase fe.g. lTornatore et al.ll2010l ). 
Unfortunately, the physical properties of both phases 
make their observ ational study challenging. Customary di- 
agnostics are Lya (|Rauchlll99g ) and metal absorption lines 
at UV and X-ray energies in th e spectra of high-redshift 
quasars (e.g. iRichter et al.l 120081 '). These tools can be fully 
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exploited only if one understands the interplay between 
galaxies and their physics, on the one side, and the bary- 
onic matter surrounding them on the other side, both in 
theory and observations. 

Ram pressure stripping, galactic winds, AGN out- 
flows and galaxy-galaxy interactions contribute to the gas 
transfer betwee n a galaxy and its ambient medium (see 
ISchindler fc Diafcrio..2008l . for a review). One important and 
widely studied facet of this problem is the chemical en- 
richm ent of the intra-cluster medium (ICM; iBorg ani et al.l 
|200g). Instead, in this work we will focus on a different but 
somewhat related point, namely the injection of turbulence 
in the gas surrounding galaxy-sized haloes. 

The discussion about turbulent motions in the gas 
at galactic scales has been addressed by several authors. 
iRauch et al.l (l200l|) studied the C iv absorption in spectra 
of double quasars at high redshift {z ^ 2), thus probing 
velocity differences (likely to be associated with turbulent 
motions) on length scales from a few parsecs to a few tens 
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of kiloparsecs. With this technique, they detect a velocity 
shear on scales larger than a few hundred parsecs, up to 
70 km s~^ at a few kiloparsecs. In the same work a comple- 
mentary approach, based on turbulence motions along the 
line of sight, derived by the width of abs orption-line profile s 
of the individual C iv systems (cf. also lRauch et al.lll996l ). 
gives a root-mean-square (henceforth rms) velocity contri- 
bution of 4.7 km s~^ on sc ales of 300 pc. 

On the numerical side, lOppenheimer fc Davg l|2009r i in 
their SPH simulations were able to fit the observed prop- 
erties of O VI absorption lines only by adding in post- 
processing sub-resolution turbulence, at the level of a few 
tens of km s~^. Turbulence injected by galactic outflows in 
the intergal actic medium (h e ncefo rth IGM) has been stud- 
ied also bv Evoli fc Ferraral Il201l|), using a merger tree as 
input for a spectral turbulence energy model. In that work, 
the average turbulent Doppler parameter (to be deflned in 
Section [3. 2 |l peaks at z « 1 to values of about 1.5 km s~^, 
with maximum of 25 km s~^. 

In this paper we will study the properties of turbulence 
in the medium around galaxy-sized haloes by means of grid- 
based hydrodynamical simulations. With respect to existing, 
similar studies in this field, we will make use of a subgrid 
scale (SGS) mod el for unresolved turb ulence. This model, 
first developed bv lSchmidt et al.l (2006) and imple mented in 
an Ad aptive Mesh Refinement (AMR) code bv Maier et al.l 
l|2009r) , has been already used for studying the WHIM prop- 
erties at cluster scales bv llapichino et al.l ((20111). The results 
of that work support the idea that turbulence is driven by 
different mechanisms in the WHIM and in the ICM, namely 
by turbulence production at shocks in the former, and by 
merger processes in the latter. At cluster scales, the level 
of turbulence in the WHIM grows with time and saturates 
aro und g = 0. 

llapichino et al.l pOllr ) and the present work explore 
very different length scales, and together with the scale also 
the involved physics is rather different. Earlier in this Intro- 
duction, we already mentioned the interaction processes be- 
tween galaxies and the surrounding medium; these processes 
are important not only for the gas transfer and the metal 
enrichment of the circum-galactic medium (hereafter CGM) , 
but also as stirring agents in the medium. Turbulence around 
galaxies is expecte d to be produced by supernova-driven 
gala ctic winds (e.g., Evoli fc Fcrrara 2011i ) and/or by merg- 
ers IjRauch et ahlbOOlr ). Furthermore, the enrichment of the 
IGM with metals at low density can be explained if feed- 
back either in the form of galactic winds or AGN feedback 
is present: these mechanisms must be invoked in order to 
provide relatively good agreement with observational prop- 
erties rela ted to quasi-st ellar object (QSO) absorption lines 
statistics (|Meiksinll2d09l '). 

In our simulations we make two (over-) simplifying as- 
sumptions: only stirring driven by structure formation, 
i.e. by mergers and shocks, is considered, and the metal evo- 
lution of the gas is not treated. These two assumptions are 
closely linked to each other: in our turbulence framework, 
any model of galactic outflow should consist not only of a 
chemical enrichment model, but also of a model for injecting 
turbulence energy. This further step, which gathers turbu- 
lence, metal enrichment and galactic processes in a single 
comprehensive model, will be subject for future work. Only 
in Section r3.3.1l we perform a preliminary study of the role of 



star formation, thermal feedback and metal enrichment on 
the production of turbulence in the CGM; however, already 
at the present stage of the project, we are able to show the 
potential of turbulence modelling in this problem, and try to 
assess what is the relative role of stirring by structure forma- 
tion. Moreover, we focus mostly on the gas with moderate 
overdensity because it pro bes the bulk of the IGM at high 
redshift (e.g. lMeiksinll2009l ) and can be studied by means of 
quasar absorption lines spectroscopy. Furthermore, the high 
redshift IGM can also be used to study the galaxy/IGM 
interplay by cross-correlation studies of galaxy properties 
and absorpti on lines along line-o f-sights to distant QSOs or 
galaxies (e.g. ISteidel et al.ll2010l ). 

This paper is structured as follows: in Section [2] we in- 
troduce the numerical techniques, the simulations and the 
gas phases that will be mainly analysed. In Section|3]the re- 
sults are presented, in terms of general properties of the runs 
(Section 13. If) , subgrid turbulence in the different gas phases 
(Section 13. 2p . in the vicinity of massive haloes (Section 13. 3|l 
and using resolved velocity diagnostics (Section I3.4|l . The 
results are finally discussed and summarised in Section O 



2 NUMERICAL TOOLS 

The simulations performed in this work were run using the 
grid-ba sed hybrid (N-body plus hydrodynamical) code enzo 
(vl.O) l|0'Shea et al.ll2005l n We simulated the evolution of 
a computational box with the comoving size of 10 Mpc h~^ 
on a side, starting from the initial redshift z = 99 (with 
the transfer functions bv lEisenstein fc Hulll999^ ■ to 2 = 2. 
Our cosmologica l parameters are, acco rding to the five-year 
WMAP results (|Komatsu et al.ll2009l ). ^a = 0.721, D,n, = 
0.279, ^b = 0.046, h = 0.7, as = 0.817, and n = 0.96. 

The computational box is resolved with a root grid of 
400"^ cells, and the same number of N-body DM particles 
(each of them with a mass of 1.01 x 10*^ M© h~^). enzo 
is an AMR code, but this tool has not been used in this 
project. This choice was made because our work is focused 
mainly on mildly overdense gas, with a baryonic overden- 
sity 5 — p/(r2bPcr) in the range between 1 and about 100 
(see below in this Section). Here p is the gas density and 
Per = SHq / (SttG) (1 + z)^ is the critical density at redshift 
z. Properly resolving this gas phase with AMR would be ex- 
tremely volume-filling and therefore very inefficient from a 
computational viewpoint, while the adopted approach per- 
forms comparatively better, with an affordable use of com- 
putational resources. The main drawback is the relatively 
coarse resolution inside the haloes. With this choice, the 
spatial resolution of our simulations is 25 kpc h~^ . In a 
similar context, the same refinement s trategy was used in 
grid-bas ed simul a tions of the WHIM bv lRegan et al.l (|2007l ) 
and Smit h et all (120111 ). 

Four simulations have been compared in our study. 
They follow the evolution of the same realisation of the ini- 
tial conditions, but use different sets of physical prescrip- 
tions. The first one is dubbed NR for "non-radiative", in 
the sense that it does not include any additional physics in 
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Table 1. Summary of the four simulations performed in the 
present work, showing the physical modules used in them. All 
runs follow the evolution of a cosmological volume with size of 
10 Mpc h~^ on a side with 400^ cells. The UV background use d 
follows the implementation describe d inlHaa rdt fc Madaul 1 11996^ , 
the cooling is treated according to iKatz et al.l 1119961). w hile for 
the turbulence SGS model we refer to lMaier et al.l ll2009h . 



Simulation SGS turbulence model 



Cooling and 
UV background 



NR 
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F 
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HC 


X 


FHC 


/ 



X 
X 

/ 



solving the equations of fluid dynamics, whose evolution is 
thus driven solely by gravity. 

The second run is labelled with F for FEARLEsfl in- 
dicating with this acronym the u se of turbulence su bgrid 
scal e (SGS) model introduced by Maier et al.l (120091) (see 
also ISchmidt etHI l2006l : llapichino et al.l I2OI0I ! I2OIII) . The 
turbulence SGS model is based on the scale decomposition 
of the fluid dynamics equa tions in a resolved and a small- 
scale part (|Germand 1 19921 ). A heuristic model is used to 
predict the kinetic energy at the unresolved, subgrid scales 
(also called SGS turbulence energy) , and to couple it to the 
resolved scales, through additional terms in the equations 
of fluid dynamics, stemming from the filtering procedure. 
For the full set of equa tions and other detail s we refer to 
iMaier et aj] (|2009l ) and llapichino et all (|201lh . One of the 
main features of fearless, namely the adjustment of the 
energy budget to the changes in spatial resolution of the 
grid, is not used in this work (with the exception of the runs 
presented in Section [33TTJ , because the computational grid 
is static and not adaptive. We also address the reader to 
iMaier et al.l (|2009l ) for numerical tests about the effect of 
resolution on the SGS turbulence, in simulations of forced 
isotropic turbulence in a periodic box. 

The third run is indicated by HC (for "Heating- 
Cooling"), because it implements th e equilibrium cooling 
model described in lKatz et al.l lll99a). togeth er with a uni- 
form U V background ( Haardt fc Madaulll996l ) . The routines 
to implement the effect of radiative cooling and UV back- 
groun d come originally from the GADGET-2 code (|Springell 
[2OO5I) and were added into e nzo for the study of the Lyman- 
a forest performed by Reg an et al.l ([2007). The He heating 
rates have been multiplied by a factor 3.3 in order to obtain 
an IGM tempera ture more in agreement with observations 
l|Viel et al.ll2009l ). No feedback mechanism is used in this 
work, but the cooling catastrophe at the centre of halo es is 
avoided by imposing a pressure floor as in Machacek e t al.| 
l|200ll ). Finally, in our fourth run, labelled FHC, the turbu- 
lence SGS model, the UV heating and the radiative cooling 
routines are used together. The summary of the runs per- 
formed is given in Table [T] 

The analysis in this work will be focused mainly on gas 
phases with selected features. Even without studying the 
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process of metal enrichment, one can infer a link to obser- 
vations by probing the gas with properties (density, tem- 
perature) which are more relevant for a comparison with 
the observables. The question at this point is: which gas 
phases are trace d by the ava il able observables? Following 
the arguments of iRauch et al.l (|200ll ). it turns out the C iv 
absorption traces gas with overdensity from a few times the 
mean density of the universe to about 200, namely a phase 
that is representative of the gas surrounding galactic haloes. 
On the other hand, the O vi absorber s traces gas with a 
slightly lower overdensity (p/p < 100; lOppenheimer et al.l 
l2012l : iTepper-Garci'a et al.ll201ll ). Without attempting a de- 
tailed comparison of our results with the above mentioned 
investigations, we will analyse mainly two gas phases: 

• a diffuse component, which in the following will be sim- 
ply referred to as IGM, with baryon overdensity 5 between 1 
and 50. We are well aware that, in literature, the term IGM 
is used in a more generic way than in our present definition; 

• a WHIM component, with temperature T between 10^ 
and 10^ K. This component is representative of shock-heated 
gas, located closer to collapsed structures. 

As we will see in the next Section, these two definitions do 
not divide the baryonic matter in an exclusive way: there 
is gas which can belong both to the IGM and the WHIM 
phase. 



3 RESULTS 

3.1 General features of the simulations 

For a first contact with the simulation data, in Figs, [l] and [2] 
we present projections at z = 2, centred on the most massive 
halo evolved in the computational box, for the simulations 
F and FHC. The runs NR and HC present only minor 
morphological differences with respect to the runs F and 
HC, respectively, thus their projections are not shown here. 

A prominent effect is the sizable role of the UV heating 
and radiative cooling in shaping the gas evolution. Simu- 
lations HC and FHC have a smaller number of filamen- 
tary structures, when compared with the NR or F simula- 
tions. This is the consequence of the efficient cooling, which 
makes gas condensing in galaxies, thus removing it from fila- 
ments. In a similar way the denser structures corresponding 
to haloes are affected, with cooling making them more com- 
pact and reaching higher densities. 

In Fig. [3] we show the gas belonging to the IGM and 
WHIM phase separately, for the run FHC. Clearly, the IGM 
gas is associated to the filamentary structure and to the dif- 
fuse medium, while the WHIM consists of gas surrounding 
more massive halos and denser filaments. The plot shows 
also that the IGM and the WHIM are not separated by 
their definitions, i.e. there is IGM gas also belonging to the 
WHIM phase. This can be seen even more easily in the mass 
distribution function of the runs HC and FHC, shown in 
Fig. 131 for all gas (black lines) and the WHIM phase (red 
lines). A significant fraction of the WHIM belongs to the 
IGM, but the peak of its mass distribution is found at rel- 
atively large overdensities {S ~ 20), consistently with what 
observed in Fig. (3) 

It has been mentioned that the use of the SGS model 
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Figure 1. The panels show a projection with size of 5 Mpc h~^ on a side, centred on the most massive halo of the computational domain 
{M200 = 1-33 X lO"*^^ Mq h~^), at the final redshift z = 2, for the simulation _F. From left to right, the panels show baryon overdensity, 
gas temperature and SGS turbulence energy. The quantities are colour-coded, according to the colour bars shown under the panels. The 
projections of density and temperature for the run NR do not differ significantly from the ones shown here. 
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Figure 2. Same of Fig. [T] for the run FHC. The projections of density and temperature for the run HC do not show significant 
differences from the ones shown here. 



does not have any significant impact on the global dynam- 
ics of the structure formation. The reason for this can be 
understood from the Mach number of the SGS turbulence, 
which is defined as 



Mt = 



(1) 



where q = (2et)^'^ is the SGS turbulence velocity, et is the 
SGS turbulence energy, and Cs is the sound speed. As we will 
see in more detail in Section 13.21 on average value of Mt is 
small (between 0.02 and 0.3) for the gas in the overdensity 
range 1 < 5 < 50. Interestingly, the average Mt is larger 
for the WHIM ((Mt) = 0.43; cf. Fig. |6l lower panel), so we 



expect that turbulence (and its modelling) are more relevant 
for this gas phase. We will come back in the next Sections 
on this key point of our analysis. 

A more quantitative comparison of the performed simu- 
lations is provided by the two-dimensional mass distribution 
function of temperature as a function of baryon overdensity, 
in Fig. El 

Here a difference between run NR and F can be seen, 
because the latter contains more mass at T < 10* K. This 
is a side effect of the SGS model, that heats unphysically 
some pristine cold gas. From equation ([T| it follows that 
in the cold gas, which has very low sound speed, even a 
moderate SGS turbulence may result in a large Mt. This 
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Figure 3. Projections of baryon overdensity for the run FHC at z = 2, lilte the left panel in Fig. [2] but showing only gas belonging to 
the IGM (left) or to the WHIM phase (right panel) For sake of comparison, the colour bars have the same range as in Figs. [T] and O 
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Figure 4. Mass distribution function for all gas (black lines) and 
the WHIM phase (red lines) as a function of baryon overdensity. 
The solid lines indicate the FHC run, and the dotted lines the 
HC run, both at z = 2. The dashed horizontal lino marks the 
overdensity range of the IGM phase. 



is indeed the case for the gas at low temperature in the 
F run, as it can be better seen from the two dimensional 
mass distribution of Alt as a function of T in Fig. [G] In the 
upper panel, a substantial amount of gas at low tempera- 
ture has Mt = 2"'^, which is the threshold value allowed by 
the safeguarding m echanism implemented in the SGS model 
IjMaier et al.l 12003 ) . which does not provide reliable results 
for larger values of Mt. The UV heating included in the run 
FHC completely solves the problem, by removing the un- 
shocked gas from the cold phase. The SGS turbulent Mach 
number still reaches very large values but only in the densest 
part of the halo cores, and we are not concerned with the 
analysis of those under-resolved regions. 

The two-dimensional T-5 mass distribution function of 
the run FHC is not included in Fig. (5] because the differ- 



ences with the run HC cannot be easily appreciated with 
this kind of diagnostic. A better visualisation is provided by 
the average temperature as a function of 5, shown in Fig. [T] 

We observe that the average temperatures agree within 
10 per cent, with the exception of the two fluctuations for 
the WHIM at 5 < 1, which have a negligible impact be- 
cause of the very low mass of the WHIM gas at that density 
(cf. Fig. |4]). In the overdensity range 5 < 5 < 40 the average 
temperature for the run FHC is slightly larger, especially for 
the IGM but also for the WHIM. A similar increase in tem- 
perature, due to the effect of including subgrid turbulence, 
has be en already observed in the simulations of iMaier et al.l 
('2003 ). and is caused by the turbulent dissipation, which 
convert a part of the unresolved kinetic energy to internal 
energy. 

The peak of the WHIM mass distribution lies in this 
overdensity range. This is probably related to the fact that, 
in Fig. [l] the WHIM mass fraction around the peak is 
slightly larger in the FHC run than in the HC. For this 
reason, it is not surprising that the SGS model is especially 
relevant for the energy budget of gas in that phase. 

At overdensities between 5 = 40 and 100 the behaviour 
is opposite, and the run HC has an average temperature up 
to 10 per cent larger. In Section [3. 3 1 this will be interpreted 
in terms of the buffer effect introduced by the turbulence 
SGS model, and examples of both this temperature trend 
and the previous one are discussed. 

Finally, the properties of the average temperature at 
even larger overdensities are more complicated. Most of the 
gas with 5 > 100 belongs to haloes, whose properties are 
not addressed in this study. An analysis will be performed 
locally in the vicinity of a few selected haloes, in Sections 
I3.3l and l3.4l Further properties of the SGS turbulence energy 
are presented in the next Section. 
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Figure 5. The three panels show the two-dimensional mass distribution function of gas temperature as a function of baryon overdensity, 
at 2 = 2. The simulations are indicated by the labels in the lower right corner of each panel. The panel for run FHC is not shown, 
because it is basically identical to that for run HC. 



3.2 Properties of subgrid turbulence 

We describe here the main features of the SGS turbulence, 
as inferred from run FHC. In this computational setup, by 
subgrid scales we mean length scales just smaller than the 
spatial resolution (25 kpc h~ ). If not otherwise specified, 
the analysis will not be extended to run F, because FHC in- 
cludes the additional physics modules (cooling, UV heating) 
which are thought to be important for a plausible descrip- 
tion of the gas physics at these scales. 

In Section TS. li the turbulence Mach number Aft was in- 
troduced for a first, brief introductory discussion of unre- 
solved turbulence. In the following we define a further turbu- 
lence diagnostic, the ratio between the SGS turbulent pres- 
sure and thermal one pt/pthcrm- Defining the SGS turbulent 
pressure as 



Pt = -pet 



the ratio is then given by 

Pt/pthorm 



2 



peint(7- 



(2) 



(3) 



which, for a value of the ratio of specific heats 7 = 5/3 
in the equation of state used at the denominator, reduces to 
the ratio between the SGS turbulent energy and the internal 
energy. Moreover, from a comparison of equations HJ and 
(O, it follows: 



Pt/Pthcrr 



7(1^^2 



(4) 



The two-dimensional mass distribution function of the 
pressure (or energy) ratio as a function of the baryon over- 
density is shown in Fig. [S] The value of the pressure ratio is 
generally low at small 5 and grows for 8 > 10. Additionally, 
two prominent features can be seen both at low and high 5. 
The first one is related to the problem of the unphysical SGS 
turbulence production in low-density gas, discussed in Sec- 
tion [STT] Here, again, the flattening of the peak is caused by 
the s afeguarding mechanism of the SGS model (Maier et al. 
l2009f ). We notice that the mass of the involved gas is safely 
negligible, in comparison with run F (Fig. [6l upper panel). 
The flat peak at higher density is related to the gas at the 
centre of halo cooling cores, and in particular to the pres- 
sure floor introduced there. This feature does not concern 



our analysis, because the innermost regions of the haloes 
will not be studied in this work. 

The average values of the pressure ratio as a function 
of baryon overdensity are plotted in Fig. [S] For the IGM 
the mass-weighted average of pt/ptherm grows from 0.002 
to 0.08. At the upper edge of the overdensity range, the 
SGS turbulent pressure starts being globally sizeable and, 
locally, there are regions where pt/pthorm can reach the value 
of 0.4 (Fig. [8|. The WHIM phase has a markedly different 
and almost flat trend with overdensity, with typical values 
of pt/pthcrm atouud 0.1. As already anticipated in Section 
3.11 turbulence is generally important in the WHIM. The 
reason for this is that low-density WHIM is contributed by 
gas that is shock heated to warm-hot temperatures in the 
outer regions of filaments and haloes, as a consequence of 
its diffuse accretion onto such structures. This is the regime 
where turbulence is injected most efficiently, thus making its 
impact on the WHIM more pronounced than in the gas at 
comparable density, but lower temperature. 

In Fig. [2] a further diagnostic of turbulent gas motions 
is introduced for its relevance in the comparison with obser- 
vations, namely the turbulent Doppler parameter bt , d erived 
from the SGS pressure like in lEvoli fc Ferrara (|201ll ): 

r- (5) 

P ^/3 

In our case, bt probes length scales at or just below the 
spatial resolution of the simulations (25 kpc h~^). The tur- 
bulent Doppler parameter 6t for the IGM grows from 0.5 
to 10 km s~^, with a volume- weighted average value of 
1.02 km s"\ while for the WHIM the average is 18.2 km s"^ 

Although at high overdensity the same average values 
are reached in Fig. |9]for IGM and WHIM, the different mass 
distributions in Fig. 2] make the large values in the IGM not 
relevant on average for this phase. 

The values of bt seen above have to be put in the proper 
context, by comparing them with the thermal broadening 
fethcrm ~ {2k-BT /m-nY i where ks is the Boltzmann con- 
stant and rriH is the mass of hydrogen, the atom that we are 
considering for simplicity. In this case, the Doppler broad- 
ening ratio can be expressed as: 



bt 



Ptmu 



Pt 



0.91 



Pt 



(6) 



fethcrm V '^'PksT ^ 'Ifl \J Ptherm V Pthcrm 

where we made use of the expression of the thermal pressure 
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Figure 6. Two-dimensional mass distribution function of the 
SGS turbulence Mach number Mt as a function of the gas tem- 
perature T, for the runs F (upper panel) and FHC (lower panel), 
at 2 = 2. 



for the ideal gas pthcrm = pTkB/(i-imp), with the approxi- 
mation for the proton mass nip ~ ttih; /i = 0.6 is the mean 
molecular weight in a.m.u.. 

From the two-dimensional distribution function of 
&t/&therm (Fig. UOP it is evident that the Doppler broad- 
ening ratio is much larger in the shock-heated WHIM gas 
than in the low-density photoionised phase. This is fur- 
ther confirmed by the averages shown in Fig. 1111 Given 
the strong link emphasised in equation ((6]), the overall 
trends for WHIM and IGM are similar to the pressure ra- 
tio in Fig. [9l The volume- weighted averages of 6t/6thorm are 
0.034 and 0.26 for the IGM and WHIM phase, respectively. 
In Fig. [TT] we show that the increase of &t/&therm with 5 
is fitted, in the IGM overdensity range, by the function 

6t/&therm = 0.023 XfJ^'^^ 

The analysis in the vicinity of massive haloes, which will 
be presented in the next Sections, concerns mostly WHIM 
gas, for which the effect of modelling unresolved turbulence 
by a SGS model are the most relevant. 
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Figure 7. Upper panel: mass-weighted average of the gas tem- 
perature as a function of the baryon overdensity 5 at z = 2, for 
the runs FHC (solid lines) and HC (dotted lines), for all gas 
(black lines) and the WHIM phase (red lines). Lower panel: ra- 
tio between the T average of the run HC and FHC for all gas 
(black line) and the WHIM (red line). The dashed horizontal line 
indicates the IGM overdensity range. 




Figure 8. Two-dimensional mass distribution function of the 
pressure ratio pt/ptherm as a function of the baryon overdensity 
5, for the run FHC, at z = 2. 



3.3 Properties of turbulence in the vicinity of 
selected haloes 

We present here a qualitative comparison of the thermal 
and turbulent properties of the gas around the second most 
massive halo formed in our computational box (chosen for 
exemplifying reasons), with mass of 1.06 x 10^^ Mq h~^ at 
z = 2, for the runs HC and FHC. This analysis provides a 
schematic picture of the most important effects of the SGS 
model on the gas physics in the outskirts of galaxy-sized 
haloes and, in prospect, of turbulence driving in this frame- 
work. 

Figs. [12] and [13] show a series of slices of gas density, 
temperature and (for the run FHC) SGS turbulence en- 
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Figure 12. Slices normal the x-axis, with size of 1.5 Mpc h~^ on a side, centred on a massive halo at 2 = 2, for the run HC. The panels 
refer to baryon overdensity and temperature, from left to right, and are colour-coded according to the colour bars on the bottom. Black 
density contours are overplotted on both panels. 
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Figure 13. Same as Fig. 1121 but for run FHC. A panel for the SGS turbulence energy is added at the right-hand side. The regions 
enclosed in the white ellipses have a temperature larger than in the run HC , as discussed in the text. 



ergy. The white ellipses in Fig. [13] put into evidence regions 
where, in comparison with the run HC, the temperature is 
higher. In two cases, the temperature increase occurs in re- 
gions where an inflowing filament (from the upper left to 
the lower right corner in the slices) drives shear flows in the 
surrounding WHIM gas. This process produces turbulence, 
as ob served in simulations at cluster scales by I Paul et al.l 
l|2011l ). For the central ellipse, a visual inspection suggests 
that turbulence is injected in that region mostly by baro- 
clinic and compress ive effec ts at the curved external shock 
(e.g.. Ilapichino &: Briiggcn |2012| V In all these cases, mod- 
elling the unresolved turbulence brings to a higher temper- 
ature in the region ar ound the inflowing filament or in the 
post-shock region. In iMaier et al.l (|2009l ). it is argued that 



the increase in temperature is caused by the additional dis- 
sipation of SGS energy. 

On the other hand, the use of the turbulence SGS model 
can also result in a temperature decrease of the gas. In 
Fig. 1131 the most apparent case is in the surroundings of 
the halo centre. This example illustrates the role of the SGS 
model as an energy buffer inserted between the resolved ki- 
netic and the internal energy component fsee llapichino et al.l 
|201(J for a graphical sketch of this concept). According to 
this interpretation, part of resolved kinetic energy is retained 
as unresolved SGS energy and is not directly dissipated to 
internal energy, resulting in some locations with colder gas 
in the simulation FHC. 

To sum up, a simple explanation would proceed as fol- 
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Figure 9. Average of the pressure ratio pt/ptherm (solid lines, y- 
axis on the left-hand side) and of turbulent Doppler broadening 
6t (dotted lines, y-axis on the right-hand side) for all gas (black 
lines) and the WHIM gas (red lines). 




Barvon Overdensity 

Figure 10. Two-dimensional distribution function of the average 
of the Doppler broadening ratio fet/fethorm iii the T — 5 plane. 



lows: it is conceivable that soon after turbulence injection, 
the inclusion of the turbulence SGS model would cause a 
reduction of temperature. On the other hand, subsequent 
dissipation of SGS turbulence could cause instead an in- 
crease of temperature with respect to the case in which no 
subgrid turbulence is included. It would be interesting to un- 
derstand under which conditions one of these two effects pre- 
vails, and to find a clear correlation between the variations 
of temperature in the runs HC and FHC and some quantity 
related to the SGS turbulence, such as the SGS energy, or 
any s ource term in its governing equation (|lapichino et al.l 
[201l|). Several tests (not included here) have not shown any 
correlation of this kind; this result is not surprising, given 
the highly nonlinear nature of the structure formation prob- 
lem. For this reason, we will not try to bring this analysis 
beyond the mere qualitative level proposed in this Section. 
Further work, based on idealised simulation setups which 



Figure 11. Averages of the Doppler broadening ratio fet/fetherm 
for all gas (black solid line) and the WHIM gas (red line) as a 
function of the baryon overdensity. The dashed line shows the fit 
in the IGM overdensity range fet/bthcrm = 0.023 X 5°-^*, slightly 
shifted downwards for case of visualisation. 



are of simpler interpretation with respect to full cosmologi- 
cal simulations, might be able to shed light on this point. 



3.3.1 Effect of thermal and chemical feedback 

As already discussed in the Introduction, the simulations 
presented in this work do not include the effect of star forma- 
tion and feedback. In this Section we will briefly report some 
results of runs in which modules for star formation, thermal 
feedback and metal enrichment are used in the simulation 
code. These results are intended to be an intermediate step 
towards a motivated implementation of kinetic feedback (at 
resolved and SGS scales) from galactic outflows, which is 
left for future work. 

Different from the runs presented elsewhere in this pa- 
per, for the simulations shown here we used the enzo code 
in its version 2.1, plus modiflcations for the implementation 
of the FEARLESS tool. The cosmological parameters, box size 
and initial/flnal redshifts are the same as in Section [5] For a 
better resolution around galaxy-sized halos, the root grid is 
resolved with 128^ cells and the same number of N-body par- 
ticles, and AMR up to three additional levels was used (effec- 
tive spatial resolution of 9.7 k pc fe~^), with refinement cri - 
teria based on overdensity (cf. Ilapichino fc Niemeveij 120081 ) 
and refinement thresholds of 8 and 4 for DM and baryonic 
matter, respectively. These choices result in a spatial reso- 
lution for the WHIM gas that is comparable or better than 
in the runs presented in Section |21 for the IGM gas, the 
same is true for 5 > 10. Like in t he previous runs, f or the 
UV background the prescription of iHaardt fc Madaul (|l996l ) 
was followed. 

Star formation and stellar thermal feedback (internal 
energy returned to the baryonic gas) was modelled accord- 
ing to the approac h and the standard parameters used by 
ISmith et al.l (120111 ). As for the metal feedback, two different 
strategies have been tested: 

• together with thermal feedback, metal enrichment to 
the CGM is also considered, with the same scheme, efficien- 
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Figure 15. Averages of the turbulent Doppler broadening fet as a 
function as the baryon overdensity (5, for the simulations described 
in Section [3.3. II (indicated by the legend). The group of black lines 
is for all gas, the red lines for the WHIM. 



cies and yields as in lSmith et al.l (|201ir ). In this case (which 
will be dubbed as run ThermMetal), the radia tive cooUng 
is provided by a primord ial chemistry network (jAbel et al.l 
ll997l : lAnninos et al. 199T) coupled to the metal cooling pre- 
scriptions of [Glover & Jappsen (2007); 

• thermal feedback is used but, for a more consistent com- 
parison with the other runs of this work, metal feedback is 
neglected. In this case, radiative c ooling is computed usin g 
the analytical approximation from lSarazin fc White! ([1987|), 
assuming a fully ionised gas with half the solar metallicity. 
This run will be indicated as Therm. 

In Fig.[l4]these two feedback runs are compared with a 
simulation, labe lled as FHC, in which only UV heating, ra- 
diative cooling (jSarazin fc WhitdllQSTl ) and the turbulence 
SGS model are employed. The main conclusion that one can 
draw from the cross-comparison is that the thermal feedback 
alone has a limited impact on the CGM. Indeed, an analysis 
like the one performed in Fig.[5]does not show any significant 
difference in the average of fet as a function of 5, between the 
runs FHC and Therm fFig. fTS)) . 

Already a simple visual inspection of the central panels 
in Fig. [TJ] suggests that the effect of metal feedback is more 
sizable. In particular, in run ThermMetal metal cooling has 
a relevant impact on the morphology of the external shock, 
which moves much closer to the halo centre. This effect is 
due to the enhanced cooling rate contributed from metal line 
cooling. Since metals are not efficiently spread in the CGM 
by thermal feedback, this enhanced cooling mainly affects 
gas in the more central regions, where metals are in fact 
confined. 

We decided not to include metal injection in our stan- 
dard runs, for the following two reasons: firstly, the treat- 
ment of metals is not complete without a model for galactic 
winds, and such a kinetic feedback model is currently not 
implemented in the ENZO code. Secondly, we see in Fig. [15] 
that, Eis a consequence of this lack of modelling of metal 
spread, the average of &t as a function of 5 does not change 
much between the runs FHC and ThermMetal, despite of 
the morphological differences. We leave the potential role of 



metals in this problem as potentially interesting, and worth 
to be taken into consideration for future work. 



3.4 Analysis of physical properties around 
massive haloes 

We now look into physical quantities related to volumes 
of 1 (Mpc h~^)'^ centred around each of the three most 
massive haloes, for all simulations. The following analy- 
sis is based entirely on the volume centred on the most 
massive halo formed in the computational box (M200 = 
1.33 X 10^^ Mq h~^ &i z — 2), results for other haloes be- 
ing very similar. Rather than analysing physical quantities 
related to subgrid scales, as done in the previous Sections, 
we extract the gas temperature, density and peculiar veloc- 
ity fields from the simulated grid (i.e. we do not interpo- 
late physical values at grid points, but use the quantities 
as already provided by the simulation grid). In this way 
we compare more quantitatively the physical state of the 
flow around these haloes in a relatively large cosmological 
volume. The analysis of the resolved velocity, in particular, 
complements the study of the SGS turbulence, and can sup- 
port the results obtained in Section [3. 21 

In Fig. [16] we show the probability distribution func- 
tions of the gas density (logarithmic units for the baryon 
overdensity), temperature and modulus of the peculiar ve- 
locity field in the left, middle and right panel, respectively, 
as extracted from the grid around the most massive halo. 
The differences between the runs that include and do not 
include UV heating and gas cooling are rather large, as we 
expect from the conclusions drawn from the T — 5 diagrams 
(Fig. [5}, while the differences between the FHC and HC 
runs are very small. 

In Fig.[l7]we consider instead the peculiar velocity fields 
and we bin the values as a function of the distance from the 
centre of mass (COM) of the halo. The velocity fields are 
particularly interesting since they are expected to impact 
on the properties of QSO absorption lines. In the j/-axes we 
show the rms value for the peculiar velocity along the x- 
direction (left panel), the rms value for the modulus of the 
peculiar velocity, in the middle panel) and the mean value 
for the modulus of the peculiar velocity in the right-most 
panel, all in km s~^. 

We decided to consider three different physical quanti- 
ties (that are of course related) in order not only to highlight 
the effects along the line-of-sight (root mean square values 
for the peculiar velocity along the a;-direction, a quantity 
that is related to the broadening of the absorption lines), 
but also to consider physical three-dimensional quantities. 

One can see that the there are only small differences be- 
tween the four models at distances larger than 0.5 Mpc h~^ . 
The F model results in a slightly larger rms values for the 
peculiar velocity along the i-axis and for the mean modu- 
lus of the velocity as compared to the NR run, but these 
differences are at the 20 per cent level. At distances below 
0.5 Mpc h~^ the FHC and HC models are very similar, and 
radically different from the runs without cooling, because of 
the larger velocity values, which can be noticed also in the 
right panel of Fig. 1171 In the next Section we will discuss pos- 
sible causes for this more vigorous fiow, found in the vicinity 
of haloes in runs using UV heating and radiative cooling. 
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Figure 16. Probability distribution function of the gas density, temperature and modulus of peculiar velocity field around the most 
massive halo for the four different runs investigated in this project: NR - non-radiative model, HC - run with UV background and 
cooling, F — SGS turbulence model and FHC - run with the UV background and turbulence SGS model (black, orange, blue and green 
curves, respectively). 



4 DISCUSSION AND CONCLUSIONS 

In this work, grid-based hydrodynamical simulations have 
been used for the study of the features of turbulence at 
2 = 2 in the medium around galaxy-sized haloes. In the four 
runs considered in this paper, we checked the effect of in- 
cluding the additional physics that is required for properly 



modelling the circum-galactic medium (UV heating back- 
ground, radiative cooling) and, for the first time in this 
kind of simulations, a SGS model for unresolved turbu- 
lence (jMaier et al.ll2009l ). Turbulence is analysed at subgrid 
length scales, namely just below the spatial resolution of 
25 kpc h-^. 
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Figure 17. Values of the IGM rins peculiar velocity along the x-axis (left panel), the rms peculiar velocity modulus (middle panel) 
and the mean of the modulus of the velocity around the most massive halo (right panel). The values are plotted as a function of the 
distance from the centre of mass of the halo in comoving coordinates. Results arc shown the four different runs investigated in this 
project, colour-coded as in Fig. 1161 



The analysis has been focused mainly on the difTuse gas 
component, with baryon overdensity in the range 1 < 5 < 50 
(IGM), and on the warm- hot phase (WHIM) with 10^ < T < 
10^ K. The gas in the IGM phase has a lower turbulence 
level (Fig. [9J than the WHIM. The trend with overdensity 
is also different: turbulence diagnostics in the WHIM are 
roughly constant over the whole overdensity range, whereas 
in the IGM there is a steady increase with overdensity. For 
the run FHC, the ratio of SGS turbulent to thermal pres- 
sure pt/pthaim grows from 0.002 to 0.08 in the IGM, and 
has a typical value of 0.1 in the WHIM. A similar trend 
has been found for a related quantity, the SGS turbulent 
to thermal Doppler parameter ratio bt/fcthorm, which has a 
volume-weighted average of 0.034 in the IGM and 0.26 in 
the WHIM, corresponding to values for 6t of 1.02 km s^^ 
in the IGM and 18.2 km s"^ in the WHIM. In both phases, 
large values can be reached locally (cf. Fig. [S]). In the IGM, 
the increase of the ratio fet/&thcrm with overdensity is fitted 
by the function fet/&therm = 0.023 x 5°'^^ 

Modelling unresolved turbulence in the WHIM by 
means of the SGS model may lead either to a heating of the 
gas, because of the dissipation of SGS turbulence to inter- 
nal energy, or to a cooling effect, because of the buffer effect 
introduced in the energy budget by the unresolved kinetic 
energy (cf. Fig. [7|). Examples of the two cases are presented 
in Section 13.31 where we show the location of turbulence 
production (post-shock regions and shear flows around fila- 
mentary inflows), together with the effect of modelling the 
unresolved turbulence. 

The mass- weighted average value of SGS turbulence en- 
ergy in the WHIM is (et) = 1.10 x lO" erg g"^ for the 
run FHC, and is larger than the value found in the run F 
(5.7x 10^^ erg g~^). Stirring is therefore more effective when 
cooling and heating from UV background are included. The 
analysis on the resolved velocity field performed in Section 
13.41 further supports this result. A possible explanation is 



that radiative cooling increases the density contrast of sub- 
clumps, which with their motion drive more vigorous turbu- 
lence in the medium. 

With respect to the IGM, the WHIM is more turbu- 
lent because is it mostly located in the vicinity of forming 
haloes, or in the inflowing filaments (Fig. [3|, and in our 
simulations turbulence is driven only by the formation of 
cosmic stru c ture t hrough merger or shocks. We notice that 
I Dave et al.l (|200ll ) already stated that the evolution of the 
WHIM is mainly driven by structure formation, with cool- 
ing and galactic feedback playing a less important role. It 
is thus reasonable to assume that the level of turbulence 
found in the simulations presented here is a lower limit for 
turbulence in the gas around galaxies, which other stirring 
sources can increase, even substantially. 

At the present stage of this work, it is hard to directly 
compare the results with other observational and theoretical 
inv estigations. The vel ocity shear of about 70 km s^^ found 
by iRauch et al.l l|2001f ) compares well, as order of magni- 
tude, with the turbulent Doppler broadening at 5 ~ 100 
in Fig. [§1 but the length scale of our analysis is somewhat 
larger. Good agreement can also be found with the level of 
sub-re solution turbulence applied by lOppenheimer fc Davg 
(|2009l ). although their analysis refers to low redshift (z = 
0.25). 

Among the limitations of the performed approach, cer- 
tainly the most relevant one is the lack of gas exchange 
mechanisms with the galaxies. First tests (Section 13.3. 1|) 
have shown the importance of metal enrichment, with re- 
spect to thermal feedback. Interestingly, from these tests we 
deduce that modelling galactic winds can affect our results 
in a twofold way, namely by stirring the gas and also by 
spreading metals in the gas surrounding haloes. 

A logical follow-up of this project is, therefore, to treat 
galactic winds, injection of metals and turbulence within 



© 2012 RAS, MNRAS OOO.ITHTil 



Turbulence in the circum- galactic medium 13 



the same theoretical framewor]^. When modelhng turbu- 
lence within this scheme, an advantage comes not only from 
the SGS model, but also from the use of a grid-based code, 
which is better suited than SPH schemes for the treatment of 
mixing processes (cf. lOppenheimer et al.ll2012r ). Of course, 
a computational strategy foc used also on better r esolving 
haloes has to be designed (cf. Humm els et al.ll2013l . who fo- 
cus on the gas properties around a single halo in AMR sim- 
ulations). It will be also interesting to check the current re- 
sults against improved versio ns of the turbulence SGS model 
l|Schmidt fc FederrathlbOllI ). more suitable for flows with a 
moderate turbulent Mach number, like the WHIM. 
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